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ABSTRACT 

The notion of a part of phase space containing desired (or allowed) states of a dynamical system is important in a wide range 
of complex systems research. It has been called the safe operating space, the viability kernel or the sunny region. In this paper 
we define the notion of survivability: Given a random initial condition, what Is the likelihood that the transient behaviour of a 
deterministic system does not leave a region of desirable states. We demonstrate the utility of this novel stability measure by 
considering models from climate science, neuronal networks and power grids. We also show that a seml-analytic lower bound 
for the survivability of linear systems allows a numerically very efficient survivability analysis in realistic models of power grids. 
Our numerical and seml-analytic work underlines that the type of stability measured by survivability is not captured by common 
asymptotic stability measures. 


The final publication is available at http : //www. nature . com/articles/srep29654 
or via the following DOT http://dx.doi.org/10.1038 /srep2 9654. 


Introduction 

In almost all dynamical systems applicable to the real world, the stability of the system’s stationary states (periodic orbits, 
chaotic attractors, etc.) is of key interest, because perturbations are never truly absent and initial data is never exactly determined. 
Nevertheless, the asymptotic stability of the system’s attractors ensures that we can still extract sensible long-term information 
from our dynamical models. 

Complementary to the notion of stability, one can analyse whether the system will remain in a desirable regime.^ This 
becomes important when a model represents a system that we have influence on, either because we engineer its fundamental 
behaviour, or because there are management options. We often want to design the dynamics, or our interventions, such as to 
more easily keep the system in such a desired state. Note that the desirable region not necessarily contains a stationary state. 

For the traditional notion of asymptotic stability against small perturbations, the key mathematical concept is the analysis of 
the linearized dynamics, in particular by means of the Lyapunov exponent or master stability function.^’ ^ 

Real-world systems typically are multistable.'”® They have more than one stable attractor,^ and thus potentially exhibit 
a wide range of different asymptotic behaviours. The key question then becomes from which initial state which attractor is 
reached, i.e., to determine the basin of attraction of an attractor. Most work so far focused on the geometry of the basin of 
attraction® of desirable attractors, e.g. by finding Lyapunov functions.^”" 

A recent idea that has been found to be useful is to study a more elementary property, i.e. not which states go to an attractor, 
but just how many. This quantity, the volume of the basin of attraction of a given attractor, can then be interpreted as the 
stability of the system in the face of a random, non-small perturbation. It quantifies the probability that the typically non-linear 
response to such a perturbation will lead the system to a different, undesirable attractor. This probability is called the basin 
stability (Sb) of an attractor.'^ This is important for a number of applications where relevant system deviations are typically not 
small, for example in neuro science, system Earth or power grids. 

One of the key appealing features of Sb is that, by studying just the volume rather than the shape of the basin of attraction, it 
becomes numerically tractable to analyse even very high-dimensional systems. It was also shown that the information revealed 
by the volume of the basin genuinely complements the information provided by the Lyapunov exponents of the system.*^ 
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There are, however, two major drawbacks when estimating Sb- On the one hand, the measure relies on identifying the 
asymptotic behaviour of a system, which might be difficult to detect, typically requires prior knowledge about the attractor’s 
nature, and is only meaningful in multistable systems. On the other hand, a Sb estimation is insensitive to undesired transient 
behaviour of the system, i.e. if the trajectory visits an undesired part of the phase space where the system would take damage 
that is not modelled explicitly. To detect this type of dangerous transients, a new, complementary measure is required. 

In this paper we introduce a new stability-related measure, the survivability S{t) of a dynamical system. This is the fraction 
of initial system states (i.e. arising from an initial large perturbation) giving rise to evolutions that stay within a desirable regime 
up to a given time t. The set of these initial conditions is called basin of survival. 

More formally, call the phase space of our system X, and a chosen desirable region C X. The finite-time basin of 
survival Xf C X^ is dehned as the set of initial conditions in X for which the entire trajectory over the interval [0,f] lies in X^. 
We choose a probability measure /r of initital conditions, reflecting our knowledge of the nature of perturbations we wish to 
study. Accordingly, the finite-time survivability is defined as 


S^it):=^i{xf). 


( 1 ) 


The total survivability then is the infinite-time limit of (t). This can naturally be decomposed into the probability that the 
initial perturbation is survived, and that the following trajectory stays save; 


SM(f) = ^0jM(X+)=S^v(t)-M(X+). (2) 

with := ]U(-nX^)/ji (X+). Now l-i{X^) does not depend on the dynamics but only on the desirable region and 

the perturbations, i.e. it is a constant for given The conditional survivability (f) captures the interplay of dynamics, 
desirable region and perturbations; it has a natural interpretation as the conditional probability of a system to survive random, 
large perturbations that do not kill it immediately. 

Assuming a uniform distribution of perturbations, the measure /r is proportional to the volume Vol. The resulting conditional 
survivability is our main object of study in what follows. We will call this finite-time survivability of a dynamical system; 


= (3) 

We are also interested in initial perturbations that only occur in a particular region of phase space. Thus, we want to study 
uniform perturbations in a subset C CX. The conditional survivability (f) can then simply be defined with respect to the 
measure Vol'^(-) =Vol(-nC)/Vol (C); 




Voi(v/nc) 

Vol(C) 


(4) 


An important example of such a conditionial survivability is the single node survivability for networked systems. There 
we condition on the phase space at a single node, thereby isolating the impact of local perturbations on the whole system. 
A mathematically precise discussion will follow in the power grid example in the results section and the supplementary 
information (SI). 

To further illustrate this definition, consider a simple example; A penguin wishing to ski down a mountain X going the 
fastest route possible in Fig. 1. The system is multistable as the penguin might end up in the goal or the valley. However, if the 
penguin goes over the cliff it will almost certainly slide the rest of the way to the goal on its back. The state of the penguin is 
not explicitly modelled by our (potential) landscape. We take this into account by declaring the parts of the cliff our penguin 
can not ski safely as an undesirable region. Further, if the penguin wishes to continue skiing, the valley might or might not be 
undesirable as well. Depending on these choices, different starting points can be in the basin of survival. If the goal is the only 
desirable attractor, the basin of survival lies in its basin of attraction, but if the valley is OK, too, this is not the case, and the 
asymptotic structure plays no role. 

As opposed to Sb or a linear(-ised) analysis based on Lyapunov exponents, the survivability is concerned not just with 
the asymptotic behaviour of the system, but depends strongly on the transient dynamics. As opposed to Sb it is applicable in 
unstable, mono-stable, or multistable, linear or non-linear systems. 
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Figure 1. Survivability cartoon: (Colour online.) A penguin can ski down the mountain starting anywhere on the slope. 
Starting at A the penguin will tumble over the cliffs, passing an undesirable state although ultimately reaching the goal. 
Starting at B the penguin will reach the goal standing on its feet. Starting even further to the right, it might end up in the valley, 
which might or might not be desirable. 


The application of the survivability concept is especially appropriate when interventions happen at the same time scale as 
the system dynamics, or when entering an undesirable region is deadly. 

A key insight is that evaluating survivability becomes amenable to Monte Carlo integration. This is due to focusing on 
the probability that the trajectory following a perturbation violates the boundary rather than trying to find the actual sets of 
phase space from which a trajectory survives. Hence, a survivability analysis, just as Sg, is applicable to very high-dimensional 
systems. In fact, the situation is more favourable than in the case of Sg, as the entire curve Sit) can be evaluated at a 
computational cost not exceeding that of Sg, while potentially revealing much more information. 

This sets survivability apart from formally similar approaches, e.g. in control theory. Their precise relationship to 
survivability is discussed in detail in the Methods section. 

For linear systems with a polyhedral desirable region, we derive a closed form lower bound on the infinite-time survivability 
5oo := 5(f —?► oo) as well as a semi-analytic, stronger bound that becomes exact in the case of vanishing dissipation. These 
bounds reveal that the survivability of linear systems depends strongly on the eigenvectors of the linear dynamics, rather than 
just the eigenvalues. The semi-analytic bound eliminates the need to simulate the system trajectory opening survivability up to 
a wide range of applications for which numerically estimating the full dynamics is not feasible. 

Results 

To demonstrate the diverse applicability of our survivability concept, we apply it to three paradigmatic model systems. A 
two-dimensional model of carbon stock dynamics, a system of integrate-and-fire neurons and a high-dimensional network 
model of the power grid. 

These systems were chosen to cover a wide range of types of systems. The carbon cycle model has one or two attractors, 
depending on the parameter regime, and some transients are deadly. The neurons are mono-stable but exhibit transient 
chaos.Finally the power grid model is high-dimensional, non-linear and multistable. However, the acceptable operating 
regime is close to a certain class of fixed points, thus the linearized behaviour near these fixed points is of great practical 
importance. 

In all three systems there are externalities which are not or cannot be modelled explicitly. Namely, the influence of dramatic 
climate changes on society, external stimuli for a network of neurons and frequency control mechanisms in the power grid. We 
will see that survivability accurately captures the interplay of externalities with the intrinsic dynamics. 

Carbon cycle model by Anderies et al. 

We begin by applying survivability to a two-dimensional carbon cycle model from climate science which has been recently 
introduced.^** This is a conceptual model with the aim to reproduce the non-linear dynamics of the carbon cycle in the Earth 
system. The boundaries of the survival region are closely related to the concept of planetary boundaries.* This system exhibits 
both the property that the undesirable states are deadly and that in some parameter regimes there is only a single stable attractor 
of the asymptotic dynamics. 
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The model equations for the atmospheric (ca), marine (Cm) and terrestrial (q) carbon stocks are given by 
Cm — dm is^a (5) 

Q =NEP{Ca,Ct) - act 

where denotes the atmosphere-ocean diffusion coefficient, j3 the carbon-solubility in sea water factor, a the human 
terrestrial carbon off-take rate and NEP (c^, Cf) the net ecosystem production, a complex non-linear relationship between the 
atmospheric and terrestrial carbon stocks (see Anderies et al.^*’ for further details). Note that the total amount of carbon is kept 
constant, leaving us with the marine (cm) and terrestrial (q) carbon stocks as independent variables. 

Part of the phase space X of the model are states with virtually no terrestrial carbon, referred to as desert states. While the 
model can recover from such states and eventually reach high terrestrial carbon states again, entering a desert state would lead 
to the collapse of human civilisation and thus, tragically, our model would no longer be valid after entering this regime. Hence, 
we define the set of desirable states as the complement of the desert states plus a safety margin m: 

= {{ca,Cm,c,) GX : c, > m} . (6) 

The safety margin should at no time, during the transient or asymptotic behaviour, be crossed. The finite-time basin of 
survival, here introduced as Xf, is then given by 




Figure 2. (a) Phase portrait of Anderies’ model (Eqn. 5, a = 0.1). (Colour online.) We choose initial terrestrial (q) and 
marine (cm) carbon stocks, the colour scale then indicates the minimum of q over the whole time evolution commencing from 
a point. An example trajectory with a long excursion to the desert state (q < m) is plotted in blue and ends at the attractor 
which is circled in yellow, the stream plot indicates the vector field of the right-hand-side (cf. Eqn. 5). The dashed black line 
indicates the value of the safety margin m — 0.1. 

(b) Bifurcations in the carbon cycle model: (Colour online.) Basin stability (Sg, blue) and finite-time survivability 
(S(t — 500), green) estimates for different values of the terrestrial human carbon off-take a. For the survivability estimation we 
assumed a safety margin m — 0.1. The shading around the curves indicates one standard error, the background colour indicates 
the different dynamical regimes. In the inset, we give survivability curves for five selected values of a, i.e. 

Ct € 10.1,0.2,0.3,0.4,0.5} from top to bottom as indicated. 


A phase plane analysis for this model is illustrated in Fig. 2a. Of special importance here are those trajectories (exemplified 
by the blue trajectory in Fig. 2a) that first cross the safety margin, i.e. are not desirable due to the very low terrestrial carbon 
stocks Q, but eventually will return to the desirable region A+. These trajectories are counted for the Sb estimation, since they 
eventually approach the attractor, but are disregarded for the survivability, since they cross the safety margin during the transient 
period. 
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By varying the human carbon off-take a in Eqn. 5, the system undergoes a bifurcation changing the number of attractors 
(around a = 0.35) as illustrated in Fig. 2b. The main picture shows the asymptotic survivability, the inset contains the 
survivability curves for different values of a. We see that the survivability drops to the asymptotic plateau at around the same 
time. Thus, if a trajectory eventually leaves the desirable regime, the time it takes until it does so is not strongly affected by a. 

The bifurcation, which is known to be a saddle-node bifurcation,^® has a drastic impact on the Sb estimation, the survivability 
only changes marginally in this interval. On the other hand, the behaviour in the interval a G [0;0.35] shows how the Sb 
estimation becomes insensitive to system changes if the multistability is lost, i.e. if there is only a single attractor (in this case 
with non-zero Cf). The crucial question whether trajectories stay in a desired regime is thus not captured by the Sb measure, but 
can be answered with the survivability concept. Note that in this case and in what follows we estimate a hnite-time survivability 
for the entire simulated time evolution of the system. Given that the asymptotic behaviour sets in earlier than the simulation 
ends, this is a good estimate for the inhnite-time survivability. 

It was argued*^ that Sb can also serve as a better early warning indicator of approaching tipping points than other measures. 
Here we see that a survivability estimation mirrors the trend in the system’s behaviour, i.e. how the set of surviving states 
depends on system parameters, while Sb remains hxed at its plateau value. Hence, survivability can serve as a complementary, 
and in some scenarios better early warning sign than Sb- 

Network of integrate-and-fire oscillators 



Figure 3. Survivability curves for networks of integrate-and-fire oscillators: (Colour online.) Finite-time survivability 
S{t\p) for given survival times f vs. the network parameter p. For each value of p we average over an ensemble of 100 network 
realisations, each with initial conditions drawn at random from the full state space. 

In the case of transient chaos*^“*^ there are long, interesting transients but potentially just a single global attractor. As 
an example, we consider a network of N integrate-and-fire neurons.^’They exhibit long-term chaotic transients, but 
asymptotically have a global periodic attractor where the neurons are in a state of phase-synchronisation. Considering the 
synchronised state as undesirable, the integrate-and-hre neurons are an example of a system in which neither asymptotic nor 
basin stability are informative. 

Modelling external stimuli as essentially randomly resetting the phases of stimulated neurons, the survivability S{t) here 
carries the interpretation of the probability that the system will not fall into a synchronised state in between stimuli, spaced 
apart at interval t. Such synchronised states model epileptic seizures and are thus undesired. 

Concretely we study the convergence from arbitrary initial conditions to periodic orbit attractors, in which several 
synchronised groups of oscillators (clusters) coexist.^® In the network every oscillator j = 1... A is connected to another 
oscillator / ^ y by a directed link with probability p. A phase variable 0j(f) S [0,1] specifies the state of each oscillator j at 
time f. The free dynamics of an oscillator j is given by 


0,(0 = 1 . ( 8 ) 

The oscillators interact on a directed graph by sending pulses when they reach the threshold <j)j = 1. After a delay time T 
this pulse induces a phase jump (indicated by differentiating the left and right limit of f as f+ and t^) in the receiving oscillator 
i: 


Ut+):=U-\u{(l)iir)-eij)) (9) 

for a potential U and coupling strength etj (For more details cf. the SI). 
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The survivability 5(f |p) for a directed network of = 16 pulse-coupled oscillators in dependence on the average connectivity 
p is illustrated in Fig. 3. For each value of p we create an ensemble of 100 network realisations. The randomly chosen initial 
phase vectors for each realisation are distributed uniformly in [0,1]^. 

All different network realisations with their associated initial conditions eventually lead to a fully synchronous state. 
However, our concept of survivability reveals the highly non-linear, non-monotonic dependence on the network connectivity p. 
While the survivability of transient dynamic states is small for networks with low and high connectivity values p, it becomes 
very large for intermediate connectivities, even for only weakly diluted networks (Fig. 3). The finite-time survivability reveals a 
new, collective time scale that is much larger than the natural period, 1, of an individual oscillator and the delay time, T, of the 
interactions. 

These long, irregular transients are the main property of interest for the system, motivating their study in.^® The dependence 
of the average lifetime of the transient chaotic trajectories on p was already studied ibidem. In this example, survivability 
reveals the same dynamical information as previous studies. Note that this is due to the specific choice of desirable region as 
the non-periodic parts of state space. Generally, there is no direct relationship between survivability and transient lengths, the 
fact that the desirable region can be chosen such that survivability reveals the quantity of interest for this system in a natural 
way speaks for its universality. 

Survivability again is a natural and informative stability measure of this system, however, this time not against perturbations, 
but against getting trapped in an undesired corner of phase space. 

Power grids 

Power grids are subject to a variety of failures and perturbations and there are numerous studies concerning asymptotic stability 
analysis, e.g.,^^’^^ and recent approaches to an Sb assessment.However, contrary to common model assumptions, the 
dynamical system does usually not evolve freely after a perturbation. If the system does not return to a stable operating state 
after a typical time span of a few seconds or if predefined thresholds are exceeded, control mechanisms that would require 
independent modelling are triggered. 

The long-term behaviour and stability of the system is thus a question for control theory rather than just dynamics. 
Conversely, the transient dynamics, and the question whether there is a temporary amplification of perturbations, is critical to 
whether the control has to be activated at all, or the system is explicitly resilient to such perturbations. Hence, the power grid is 
an example where the undesirable region is deadly and management options operate at the system dynamics time scale. 

The effective network model of the power grid^*’^^ is the current standard baseline model for the frequency dynamics of 
power grids. It is known as the swing equation or the second-order Kuramoto model, and is used for short-term frequency 
stability studies in power grids. The various ways in which a power grid can be modelled using the swing equation are discussed 
in^^ and limits to its applicability are discussed, for example, in.^^’^"* 

The dynamical system modelling N generators’ instantaneous phases and frequency deviations (Ot from the grid’s rated 
frequency is given as 

= (Oi (10) 

N 

d)i = Pi - am - ^ Kij sin {(pi - <pj) 
j=i 

with Pi being the net input power/consumption, a, the electro-mechanical damping at node i and Kij as the capacity of the link 
i - j. Here we choose F) = 1 for net generators, P, = — 1 for net consumers, and a uniform distribution of a, = a = 0.1. We 
choose the nonzero Kij uniformly equal to 6, corresponding to an average transmission line length of about 200 km. 

A stable operating state of the power grid is a fixed point of the dynamics with no frequency deviation, {(p*,0) : = 
{(pl ,... ,0,...). Conversely, limit cycle solutions (frequency oscillations) need to be prevented in order to avoid the tripping of 
generators. Frequency deviations are usually kept very small in large real power grids, with typical thresholds of ±0.2 Hz^^ 
which corresponds to a phase velocity deviation of |a)| « 0.25 in our units. Smaller island grids have considerably larger 
fluctuations. As an illustrative extreme case we will consider up to 20 times larger fluctuations. For Sb assessments, the reaction 
of the system to much larger deviations was also taken into account. 

We will study the single-node basin of survival, i.e., the conditional basin of survival in the sense of Eqn. 20, conditioned on 
initial perturbations that occur locally at a single node n, starting from a stable operating state. The space we wish to condition 
on is then the direct product of the stable operating state at all nodes except node n and the full state space of the node dynamics 
at n: 


Cn — {(•/•! )■■■•/*«)•■■) 1 0,. . ., ,. 


..,0)Is 


[ 0 , 27 r),®„ e K} . 
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The desirable region being defined as Vi: |a)j | < 5, which, as explained above, is chosen to mirror realistic constraints. 
Concretely, this means that we construct initial conditions by setting 0, and m, to the value of the fixed point 0* and 0, for all 
nodes other than the node n we are smdying, and to a random phase in [—n\ n] as well as a random frequency deviation in 
[—5; 5] for the node n. Then we simulate the system up to f = 100 and observe whether (and if, when) any of the frequency 
deviations O), leave the desirable region. In this way we sample 300 trajectories to estimate 5”(f) := 5*'" (f). This leads to a 
standard error of less than 0.03 for 5"(f) = 0.5 in the worst case (see Methods section). We evaluate the survivability up to 100 
in simulation time (18s in real time), at which point a steady state has typically been established, and the asymptotic value of 
the survivability is reached. 

While Sb captures the overall ability of the system to avoid permanent frequency oscillations, it does not directly capture 
the stability of the system against large perturbations. Instead, as discussed above, it is the ability of the system to keep 
perturbations under fixed frequency thresholds which is crucial. We will study this form of stability using both numerical 
simulations and the analytic approximations we have derived. The former will allow us to compare the survivability of the 
system to its Sb, the latter to assess the accuracy of our bounds. 

We now turn to the question whether the semi-analytic bounds on the dynamics linearized around the fixed point can 
accurately mirror the single-node survivability S"(t = 100). 

Defining 0 := {(pi,..., (PnY , (0= (cui,..., COnY ^nd a := diag{ai), the linearized dynamics is given by 


/0\^/O In\ f(p\ 
J \L -aJ \(oJ 


( 11 ) 


where the lower left block (L = dati/dpj) can be identified with the network’s Laplacian matrix (at the fixed point {(p*,Q)) 
given by 


Lij = -Sij £ Ki„cos iP* - P:,)+Kijcos {P* - p*) 


m=\ 


( 12 ) 


The Jacobian has two real eigenvalues, = 0 and A2 = — Ct, corresponding to the eigenvectors {p,co)i = (1,... ,0,.. .) 
and {p,co )2 = (—l/ct,..., 1,...). The first eigenvalue, Xi and the corresponding eigenvector show the linearized version of 
the rotational symmetry of the system under shifting all elements of p by the same amount ps'. pi 1 —pi + ps- The second 
corresponds to a homogeneous shift of all oscillator’s frequencies, which does not affect the phase differences, and decays 
exponentially due to the damping term. The remaining part of the spectrum consists of — 1 pairs of complex conjugated 
eigenvalues. 



Figure 4. Single-node phase space of a consumer in the Scandinavian grid: (Colour online.) (a) We plot the initial 
frequency deviation m,, vs. the phase difference to the fixed point at node n, visualising the definition of the following areas 
using the simulation results from (b). The central green area resembles the infinite-time basin of survival, while the yellow and 
red areas contain finite-time surviving states. The union of the blue, yellow and green regions resembles the synchronous state’s 
basin of attraction, while trajectories starting in the white or red regions approach different attractors. The frequency threshold 
is chosen as COcrit. = ±5 and initial conditions correspond to perturbations at a single consumer node of the network, (b) 
Simulated maximum frequency deviations (Omax along all dimensions, measured over the time evolution of the system for initial 
conditions that correspond to perturbations at node n of the network. For comparison with (a), we give the numerically 
estimated basin of attraction’s boundaries in red. (c) Corresponding analytic upper bound for the maximum frequency deviation 
(cf Eqn. 16) for the linear approximation. 


The basin of attraction in the conditional subspace C„ of this system is illustrated in Fig. 4 (a). Concerning survivability, 
there is a subdivision in three different sets. The desirable region contains infinite- (central green region) and finite-time 
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surviving states (yellow and red regions in the band). Trajectories commencing from the remaining states within the basin of 
attraction (blue region) eventually reach the attractor asymptotically. Note that there are also finite-time surviving states outside 
the basin of attraction (red region). A large part of the single-node basin of attraction is centred around the hxed point ((j)*,0). 
Within this region we expect the linear approximation to provide a lot of information on the system. 

Regarding survivability, Fig. 4 (b) shows that the frequency deviations inside the basin of attraction do indeed become large. 
The shape of the level lines of the frequency deviations corresponds to the basins of survival for different frequency constraints. 

Fig. 4 (c) shows the bound for the frequency deviation of the linearized dynamics calculated according to Eqn. 16. This 
shows a good qualitative agreement with the actually simulated frequency deviations as long as the deviations remain close to 
the fixed point, e.g. in the range of realistically allowed perturbations (see above). Still, the impact of the non-linearity (e.g. 
multistability is not captured) on the system becomes apparent, especially further away from the hxed point. 
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Figure 5. (a) Simulated vs. approximated single-node survivability for the Scandinavian grid: (Colour online.) Scatter 
plot of the simulated S"{t) vs. approximated single-node survivability (cf. Eqn. 16) estimated for all nodes in the 

Scandinavian power grid (cOcrit. is indicated in the legend). The corresponding distributions are given on the sides. 

(b) Single-node basin stabiUty vs. single-node survivabiUty for the Scandinavian grid: (Colour online.) Scatter plot of 
the single-node basin stability Sg vs. single-node survivability 5”(f = 100) (cOcrit. is indicated in the legend) estimated for all 
nodes in the Scandinavian power grid. The corresponding distributions are given on the sides. Note that we have chosen the 
initial region Xq for single-node basin stability with |cu| < 100, the same region as in.^® 


Indeed Fig. 5a shows that there is a high correlation between the lower bound of the survivability of the linear system sUO 
calculated according to Eqn. 16 (see Methods section) and the actual survivability 5"(f) at the majority of nodes for realistic 
values of frequency deviations. What exactly gives rise to the outliers far below the diagonal will require further study. It is 
important to emphasise that the computational cost of calculating the bounds on the maximum frequency deviation for a sample 
of initial conditions is many orders of magnitude lower than the numerical estimate of the survivability via simulations of the 
actual time evolution. For a realistic network size of several hundred nodes, the approximate calculations can be performed on 
a laptop computer in less than a minute, whereas the numerical survivability estimation took several hours on 200 nodes of a 
computing cluster. 

Fig. 5b shows Sg as well as the single-node survivability of nodes in the Scandinavian power grid. We see that there is no 
signihcant correlation between the two quantities. This proves the point that the asymptotic behaviour of the system is not a 
strong indicator of the transient behaviour, at least in the case of power grids. The information we obtain from the survivability 
analysis is genuinely new information. 

The Scandinavian power grid^® consists ofN = 236 nodes and 320 links, corresponding to a mean degree ofk = 2.7. Hence, 
it has a sparse network topology with only a few neighbours per node on average, which is typical for power grids in general, 
independent from the number of nodes.The same holds for our second data set, the UK high-voltage transmission grid, 
which consists of N = 120 nodes and 165 links, corresponding to a mean degree ofk = 2.8. 
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Figure 6. (a) Scandinavian power grid: (Colour online.) The nodes’ colouring indicates the respective single-node 
survivability estimate S"{t = Is) in the Scandinavian power grid. The frequency threshold is chosen as (Ocrit. = ±10. We 
randomly selected a dispatch scenario, circular nodes are net generators, squares are net consumers. The map of Scandinavia 
has been modified from https : //commons . wikimedia . org/wiki/File : Scandinavia . svg, which is licensed 
under the Attribution-Share-Alike 3.0 Unported license. The license terms can be found on the following link: 
https://creativecommons.org/licenses/by-sa/ 3.0/. 

(b) UK power grid; (Colour online.) Single-node survivability estimate 5"(f = li) of the UK power grid. Details analogous to 
(a). The map of Great Britain has been modified from https :/ / commons . wikimedia .org/wiki/File: 

England, _Scotland_and_Wales_within_the_UK_and_Europe . svg, which is licensed under the 
Attribution-Share-Alike 3.0 Unported license. The license terms can be found on the following link; 

https://creativecommons.org/licenses/by-sa/ 3.0/. 


In Fig. 6a and Fig. 6b we show the geographically embedded Scandinavian and UK power grid. The colour of each node 
corresponds to the single-node conditional survivability 5"(f = Is). Different nodes exhibit starkly different survivability to 
perturbations. We find that at a threshold of I = 10, for both of these realistic power grid topologies, there are a few 
nodes that are particularly vulnerable to perturbations. This means a perturbation at these nodes is very likely to be amplified 
temporarily by the overall grid dynamics. What exactly leads to this vulnerability, and how to characterise it in terms of grid 
parameters and topology is a question for future work. 

Finally, we also found that the survivability in this system asymptotes very quickly. Simulating just the first second of the 
power grid is typically sufficient, the so-called “first swing” following a disturbance mainly determines the overall frequency 
deviation. 

Let us summarise the key points from applying survivability to power grids: 

• For realistic small deviations, the upper bound applied to the linear approximation provides an excellent picture of the 
infinite-time basin of survival. The fact that the bulk of nodes shows a high correlation at large perturbations indicates 
that 5" can still be determined from the approximation in this case. 

• For the given dynamics, the survivability very quickly reaches its asymptotic value. We expect this to be a fairly generic 
phenomenon if we are dealing with damped systems near a stable fixed point. 

• Conditioning the survivability on regions of phase space with special meaning, like perturbations at a single node, allows 
us to reveal a large amount of non-obvious structural information on a networked system. Further work is needed to 
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understand what gives rise to the revealed structure in realistic power grids. 


Discussion 

Survivability is a novel stability concept complementary to basin stability Sb and linear methods of asymptotic stability analysis. 
It applies to linear and non-linear systems, in the absence and presence of multi-stability. It focuses on transient rather than 
asymptotic behaviour, and incorporates exogenous information via assuming a desirable region for the system dynamics. 
Further, survivability can be estimated numerically at low computational costs, comparable to or even lower than for estimating 
Sb- 

For linear systems we provide easy to evaluate analytic and semi-analytic expressions for lower bounds of the survivability, 
with a trade-off between the quality of the bound and numerical cost for evaluating the analytic expression. These reduce the 
need to simulate the system, yielding further dramatic improvements in computational cost. 

The bounds we hnd demonstrate that the survivability depends crucially on the eigenvectors of the linear dynamics, rather 
than the eigenvalues (see discussion in the Methods section). It is an effective measure of the interaction between external 
constraints and the geometry of the dynamics in its phase space. The fact that the bound is tight exactly when the analysis of 
asymptotic stability using the eigenvalues of the linearized system fails shows that the survivability is genuinely complementary 
to eigenvalue-based stability concepts. 

To explore this measure in practice, we analyse three conceptual examples. 

Carbon Cycle: We observe that survivability accurately exhibits the presence of dangerous transient behaviour in the 
model, something that Sb can not detect. The almost monotonous decrease towards the first tipping point, opposed to the 
discontinuous Sb curve, shows the potential to derive an early warning scheme from an observation of these measures for 
certain kinds of bifurcations. Just as for Sb, the problem of evaluating the survivability from data remains a challenge for future 
work. 

Neuronal Networks: Here, the transients do not arise from perturbations constructed as deviations around a desirable 
attractor, but they are randomly chosen from the whole compact phase space. Rather, the main interest lies on the transients 
themselves. Survivability reveals the same qualitative dependence of the dynamical behaviour on the underlying network 
topology as the average length of the transient.^® Beyond that, considering S{t) at fixed f as a function of the underlying 
topological parameters enables us to look in more detail into the relationship between function and structure of pulse-coupled 
oscillator networks. In contrast to the average length of the transients, the survivability also has a direct conceptual interpretation 
as the probability of the system remaining in the interesting transient regime. Thus it captures the appropriate notion of stability 
of transient chaos against the global attractor. 

Power Grid: In this example we can see in detail the interplay between the semi-analytic bounds that we developed and 
the fully non-linear system. We demonstrate that survivability under realistic constraints captures information about the system 
not contained in the Sb estimate. We also demonstrate that the semi-analytic lower bounds, are strongly correlated with the 
simulations of the non-linear dynamics. Thus they contain much of the relevant information about the system. In strategic 
power grid development studies, this fact becomes particularly important as computational power is often at a considerable 
constraint, due to the need to simulate a wide range of divergent future scenarios of the energy transition. Dynamical properties 
outside of quasi-stationary calculations can only be taken into account if efficient estimators exist, since it is not feasible to run 
simulations. Thus our lower bounds, which eliminate the need for such simulations, potentially enable a more systematic way 
to investigate the impacts of the energy transition. In particular, the influence of changing topologies and different distributions 
of dynamical parameters on the dynamics of the power grid become computationally accessible. For the application to power 
grids, there are many more operational conditions on the system’s behaviour that we do not consider here. While not all of them 
are as amenable to analytic considerations as the frequency deviation, we anticipate that it will still be possible to find cheap 
analytic boundaries for them. The reason that we could calculate the lower bounds so easily is that the phase space geometry 
is encoded in an efficient way in the eigenvectors. This aspect will carry over to many other, more complicated exogenous 
boundaries. 

We thus have seen that the notion of survivability is general and powerful enough to capture the interplay between 
externalities and the intrinsic dynamics in three vastly different examples. In particular the last example demonstrated both the 
utility of single node survivability, revealing structural weaknesses and strengths of realistic power grid topologies, as well as of 
our semi-analytic bounds, reducing computational efforts dramatically. 

The work presented here thus opens up a plethora of new avenues of research. On the theoretical side, the existence of 
a closed form lower bound on the survivability of a linear system opens the door to study the survivability as a function of 
the network topology and system parameters analytically, especially for the optimisation of these parameters to increase the 
system’s survivability. The lower bounds presented here can certainly be improved by taking the more detailed geometry of the 
trajectories of the linear system into account. It will also be important to extend them to the types of bounds we have in more 
realistic power grid models. 
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Methods 


Numerically estimating survivability 

One advantage shared by survivability and^^ is that they can be efficiently estimated by randomly sampling starting conditions. 
A trajectory either survives or not, therefore we can regard the sampling as a Bernoulli experiment with probability given by 
S{t), hence the standard error (SE) of the probability estimator of a trial with N draws is simply 

,.3, 

As a crucial consequence, the standard error of a survivability estimation does not depend on the dimensionality of the system. 
Further, the condition that a trajectory has left tends to be easier to evaluate in practice than whether the trajectory is 
asymptotically approaching a fixed point. Furthermore, in numerical simulations, an integration might be stopped once has 
been left. 

Analytic results for linear systems 

An important analytically tractable case is the total survivability 5oo for a linear dynamic inX = the Lebesgue measure 
Vol(X) = dx^, and a polyhedral desirable region given by m linear conditions yt'xit) < 1 for a set of vectors yk, k— I ...m 
in In this case we can give a lower bound on Vol(x£) that is easy to evaluate. 

In this section we briefly give the results necessary for the applications in the results section on power grids. There we 
demonstrate that the semi-analytic bound captures the survivability of the system quite accurately in practical examples. In the 
SI we show detailed derivations, as well as further analytic results. 

Consider a system of linear ordinary differential equations 

x{t)=Lx{t) (14) 

with x&X = and L S with all eigenvalues having non-positive real parts. In general, L has a complex spectrum. The 

eigenvectors Vj of the complex eigenvalues are real or come in complex conjugate pairs, from which we pick one eigenvector 
each. We then define the N xN matrix V by stacking the eigenvectors, or their real and imaginary parts respectively, against 
each other as column vectors: 

V= [vi,...,Re(v;),...,-lm(vj),...]. (15) 

This allows us to translate initial conditions into the eigenvector basis by setting c' = V^*x(0), and combining into 
complex numbers as appropriate cj = cj, -f where tic is the number of complex eigenvalues. Then the trajectory describes 

an exponential decay along the real eigenvectors and an inward spiral in the Re (vj), Im (v;) plane that is parametrised by cj, 
and given by Re {exp{Xjt)cjVj). We then obtain an upper bound for the deviation of the trajectory starting at x(0) in a direction 
yk by maximising the contribution of each eigenvector separately. 

Now, setting ykj ■=yk‘ vj for Vj real, and ykj '■= \yk' Vj\ for Vj complex, this leads to the estimate: 


«0 fh n 

max \yk-x{t)\<yykjCj+ T max{0,ykjcj) + T ykj\cj\ (16) 

;=1 j=no+l J=nr+l 

where the first sum is over real eigenvectors corresponding to null eigenvalues, the second is over nonzero real eigenvectors 
and the last is over the complex eigenvectors. 

Setting the right hand side of Eqn. 16 smaller than 1 defines a region Vj, in spanned by the real and imaginary parts of 
the coefficients cj. This region is mapped to the state space by V and thus its volume is related to the corresponding region in 
phase space by a determinant factor. As it is defined by a weaker inequality than X^ it follows that 

Vol(Xi) > VdetW^Vol(14) ■ (17) 

The inequalities Eqn. 16 together with the matrix V can be used to efficiently estimate the total survivability as well as the 
conditional survivability. Remarkably, for systems with a purely imaginary spectmm, the bounds of Eqn. 16 and Eqn. 17 hold 
with equality. 

In the SI we also derive a lower bound for Vol(Vc). 

This lower bound demonstrates that for the survivability of a linear system, the eigenvectors play a crucial role. In fact, the 
eigenvalues do not enter the bound at all, except in terms of classifying the corresponding eigenvectors in separate classes. This 
demonstrates that the survivability captures substantially different information about the linear system than eigenvalue-based 
stability measures like relaxation time, or the master stability function. 
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Relationship to Similar Concepts 

Survivability is related to a number of concepts in other fields, notably control theory. From this perspective it can be seen as a 
so far unstudied, simplifying case where a number of distinct concepts from various fields intersect. In this section we discuss a 
number of such concepts and their precise relationship to survivability. 

Survivability is conceptually similar to the notion of finite time stability as studied for linear control systems. There the 
focus is on finding a particular control scheme that will ensure that the resulting closed loop system stays within a particular 
region for some time, possibly in the presence of perturbations of the dynamical equations. From our perspective this can be 
seen as attempting to find systems with S{t) = 1. As the focus there is on perturbed dynamics in linear control systems, the 
actual overlap of methods is very small, in particular it is not possible to extend the methods to high-dimensional non-linear 
systems. 

Another concept from control theory which is similar to the basin of survival is the viability kernel defined by Aubin et. al. 
in the context of viability theory.They introduce the notion of an environment K that contains all desirable states. Within 
the environment, there is the so-called viability kernel as the set of all initial conditions from which the system can stay 

within the environment. This basically is a more general version of our infinite-time basin of survival for non-deterministic 
systems or systems with multiple evolution paths and a management process. Consequently, K \V corresponds to the set of 
finite-time surviving states in deterministic systems. The viability kernel’s volume is proposed as a measure of the degree of 
viability,^^ in the limit of no control it thus reduces to our total survivability. However, we are not aware of this special case 
ever being considered in the context of viability theory. Whereas survivability measures the ability of the intrinsic dynamics 
to withstand perturbation, viability theory is concerned with the question of the power of control. Beyond this conceptual 
difference, evaluating survivability also requires very different technical methods, analytically as well as numerically. As far as 
we are aware, sampling based methods, which are efficient and natural for survivability, are impossible for viability. This is due 
to the fact that whether a particular point belongs to the viable set depends on the optimal control, which might not be known. 

There are two concepts that share some formal similarity to survivability in the context of deterministic systems, transient 
times and open systems. 

The study of transient life times^^’ '^’41,42 related to the survivability in the non-typical special case that the attractor 

(or a small epsilon environment around it) is the only undesirable region. In our example of integrate-and-fire neurons this is 
the case, but in the power grid there is no clear relationship between the strength of the transient (which might kill the system) 
and the return time to the attractor. In fact, there, the attractor we start from is in the desirable region. Transients life times are a 
special case, and not a typical one, of survivability. The latter is far more general, going beyond the focus on the length of 
transients and their distribution, and typically captures genuinely different information of the system (e.g. the linear analysis 
mainly depends on eigenvectors, not eigenvalues). 

The theory of open systems, on the other hand, is generally concerned with ergodic systems. For leaky chaotic systems^^ the 
asymptotic behaviour of the survival probability is the key observable. At the formal level there is an analogy to our definitions, 
however, the total survivability, the size of the total phase space that leaks, is never considered as an observable in the literature. 
Indeed it is often the case that it is the whole phase space. Nor is the cumulative leakage ever interpreted as a stability measure 
or are efficient methods to estimate it for high-dimensional systems being discussed. In fact, as in the case of transient times, 
leaky systems can be seen as a special case of our discussion. Specifically it is the conditional survivability with the conditional 
space chosen as the space of surviving states X^. 

The closest analogy to our deterministic survivability is simply the survival analysis in the the context of stochastic systems. 
The concept of the so-called first hitting time and survival probability,^^® which can be studied for the case of stochastic 
perturbations to deterministic systems by quasi-potentials,'^^^^ map directly to our work. The first hitting time t measures 
when a system is expected to first hit the forbidden region X^. The cumulative of the probability of first hitting the undesirable 
region before t is then 1 — 5(f). Our definitions given above can be seen as a deterministic version of these concepts. The role 
of stochasticity in the evolution is replaced by a probabilistic initial perturbation. Here similar sampling based methods are 
possible and necessary. The type of semi-analytic analysis we performed for the linear case would however be hard to duplicate. 
From this perspective what we have demonstrated is how to successfully apply methods and concepts from stochastic systems 
in the study of their deterministic counterparts. 

The key insight in our work, as it is for Sb, is that restricting ourselves to probabilistic notions enables a considerably wider 
applicability of our analysis, as well as new numerical and analytic methods. Put differently, by asking not about the geometry 
of sets in phase space but merely about their volume, we can access high-dimensional non-linear systems that are out of reach 
for detailed geometric analysis. The challenge then lies in defining interesting sets that capture concepts of interest. As such we 
take it as a confirmation for the wide interest of the specific sets that survivability is based on, that it occurs a the intersection of 
a number of well studied concepts. 
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A Formal Definition and Derivation of Analytic Bounds 

As noted in the main text, the survivability of a linear system is amenable to analytic study. In this appendix we will give a 
more mathematically precise dehnition of survivability and then a detailed derivation of the results used in the main body of the 
text as well as some closely related ones. 

A.1 Formal definition and basic properties 

Consider a dynamical system with states x in a state space X giving rise to trajectories x{t) under some evolution map (T(f). 
Now we dehne a desirable region A+ C X with its complement X~ = A \ A+, where the former contains all states x that are 
stated to be desirable. In the penguin example (Fig. 1) X^ would contain the cliff and the valley. In the context of Earth System 
science, such a desirable region has variously been called the safe operating space within planetary boundaries^ or the sunny 
region? 

According to the main text, we dehne the hnite-time survivability S{t) of the dynamical system at time t to be the fraction 
of trajectories starting in A+ that stay within A+ for the entire duration [0,f]. Put another way, if entering the region X^ 
terminates the system, S{t) is the fraction of trajectories starting in A+ still alive after time t. We call the part of A+ from which 
trajectories start that stay alive at least for time t the finite-time or t-time basin of survival Xf. We then have 


5(f):=V(0 = ^^, (18) 

where p is an inner measure on X determining the volume of the sets Xf and A+ in the phase space. By construction, 
(f) takes values on the unit interval. 

We dehne the total survivability (f — °o) as the limit 


5oo:= lim5(f). (19) 

Each f-time basin of survival is a subset of the previous ones X,^ D X^ (for t' > f), as trajectories returning to A+ after 
leaving it once do not contribute to Xf,. Hence, S{t) is monotonically decreasing and bounded by 0 from below, therefore the 
limit in Eqn. 19 exists. The use of an inner measure here avoids subtleties involving non-measurable sets, like fractal^’ ^ or 
riddled^’ ^ basins of attraction. 

Note, however, that if A+ is an open set, and the map ff (f) : A —A is continuous for all f, then the images of A+ under 
<j{t)^^ are also open. As Xf = Uo<i'<r <7(f')^'A+ is a union of open sets, it is itself open and therefore measurable if /r is a 
Borel measure. Thus in this important special case, which covers all applications we are considering in the main text, no such 
subtleties exist. 

In some applications, the choice for the set A+ might have an inhnite volume, even though A/ becomes hnite for sufficiently 
large t. In that case one can still consider the unnormalised measure /X (A,‘^) as a relative measure of the survivability of a system. 
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A.2 Conditional survivability 

The conditional survivability S^{t) measures the response of the system to restricted perturbations. For example, we might be 
interested in the survivability, given perturbations that are localised at a node in a network. Given a subset of the state space 
C C X, we define the conditional survivability as the fraction of trajectories starting in n C that stay in X^. That is 


^\c{x+nc) 


( 20 ) 


where /i|c is an inner measure on the smallest sub-manifold containing C. In the case of a network and perturbations at a single 
node, the phase space typically is the product of phase spaces at the nodes, and the volume measure likewise factorises, thus 
there are natural choices for C and /r|c- 


A.3 Linear Systems 

We consider the case of a linear dynamic in X = the standard Lebesgue measure Vol(X) = dx^ and a polyhedral sunny 
region given by m linear conditions yi^-x{t) < 1 for a set of vectors yi^, k = 1... m in 
The dynamics is given by a linear system of ordinary differential equations 


x{t)=Lx{t) (21) 

with x&X = and L G In general, L has a complex spectrum, and we assume that all eigenvalues have non-positive 

real part. We denote the number of real eigenvalues of L as tir, the number of pairs of complex conjugate complex eigenvalues 
ric- Then we have a total of n — rir + ric independent eigenvalues and rir -f 2«e = N. We denote the eigenvalues as A, and A,-, and 
the corresponding eigenvectors as v, and vl respectively. 

Assuming that the real and imaginary parts of the eigenvectors of L span the entire space X, the general solution to the 
dynamical equations are then given by the matrix exponential 


x{t) = expLfx(O) 

n 

= L Re(c/expA/v, ) 
j=i 

rtf n 

= ^ Cy-expAy-fV;-f ^ Re(cyexpA/fVy), (22) 

j=i ;=«r+i 

where the coefficients cj are real for j < and complex above. To determine them we introduce a convenience map i as 
follows: 


{ Cj ifj<nr 

Re (cj) if ny<j< tic (23) 

This is a real-linear map. Then we can define the real matrix V as: 


Re(t'„,+i),..-,Re(v„), 

-lm(v„,.+i),...,-lm(v„)]. (24) 

and obtain 


x(0) = (Voi)(c), (25) 

which can be readily inverted as we assumed V to have full rank. 
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A.4 Upper bound on the deviation of a singie trajectory. 

The inner product of x{t) with a boundary vector y G is then simply given by: 


rir fir+flc 

y ■ x(t) — ^ Cj expXjty ■Vj+ ^ Re {cj exp A/y • vj) (26) 

7=1 i=nr+l 

Now, to obtain an upper bound on this for all times we can maximise each individual contribution. As the eigenvalues have 
non-positive real part each complex contribution has magnitude of at most \cjy • V; |. The maximum of the real contribution 
depends on the sign of cjy- Vj and is given by 


max(0,cjyvj) (27) 

unless Aj = 0 in which case the contribution is exactly equal cjy ■ Vj. While this does not happen generically it occurs due to 
symmetries in the system, and thus we will treat it separately. We denote the number of zero eigenvalues by hq. We have: 


no 

max \y-x{t)\ < V c,y-v; 
r6[0H'^ ' " ' 

Hr 

+ E max(0,cyy vj) 

7=no + l 

+ L (28) 

7=«r+l 

This is the key approximation for our analysis. In the next section we will use this estimate to give a lower bound for the 
total survivability, afterwards we will show when this bound becomes tight. 

A.5 A lower bound for the total survivability 

For a boundary vector yj^ let us define yj^j = y^ ■ vj for j < tir and yj^j = |yj. • Vj \ for iir < j < n. Then the inequalities 

no n,- n 

'LykjCJ+ L max(0,yA:;C;) + E ykj\cj\<l (29) 

;=1 J=no + l j=nr+l 

define a region Vc in K”'' (g) C”'" that is mapped to a subset of by the linear transformation V o i. This is a subset of as the 
initial conditions in this region have an inner product with the boundary vectors y^ bounded by Eqn. 28. Thus, taking the effect 
of the transformation V into account, we have the lower bound 

Vol(xi) > \/detW^Vol(yc) . (30) 

This bound can be evaluated numerically quite easily. All that is needed is to sample A+ and take the fraction of samples 
for which the image under the linear map (Vo i)^* satisfies Eqn. 29. This also makes it very easy to numerically estimate the 
lower bound of conditional survivability, by sampling from Cn instead. The data for the semi-analytic figures in the power 
grid section of the results was computed in this way. 

In order to proceed with the analytic calculations we have to consider further special cases. First we will show that the 
bound Eqn. 30 is actually exact for some cases. 

A.6 The case of vanishing real parts. 

Let us now consider the case where the real part of all A, is zero, and 
trajectory with initial conditions x(0) = V o i o c is dense on the torus: 




is irrational. Thus no = n,-. In that case the 


u'=l \j=nr+l J 


0;e[O;27r[ 


(31) 
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The maximum of the torus along the y direction is obtained exactly by maximising each contribution independently, which 
means tuning the 0/ so that Cjexpi(j)jX ■ vj are real, thus the real part equals the absolute value; 


max \yT\<Y^Cjyvj+ £ \cjy 


c|i,e[0;2;r[ 


(32) 


;=i 


;'=«(■+1 


with equality if is irrational. Therefore, in this case, the general bound Eqn. 28 holds with equality. We have; 

I m ( Ay j 


rir n 

max |yx(f)| = £c,-(vo)yV;+ £ |c;(xo)||y V;| (33) 

;=i j=nr+l 

and therefore also 

Vol(Xi) = VdetW^Vol(W) . (34) 

A.7 The purely imaginary case. 

For the case = no = 0 we can give an explicit lower bound for Vol(yc). To do so, define yj = miny^^j. The volume of the 

k 

space Vc defined by the inequality 


E>vk./I < 1 

./=i 


(35) 


is a lower bound for the volume of V) . The two spaces have the same volume if there is one condition that dominates all 
others, that is, if there is a k' such that yj = for all j. This means that any trajectory that leaves the allowed region also 
crosses the boundary defined by x y^; < 1. 

We now need to evaluate the 2n = 2nc = N dimensional integral 


Vol(l4) = / Y\Ac)dc) , 




with Cj = c'” + ic‘j. We begin by changing dcydc' to polar coordinates rjdrjdipj and rescaling; 


(36) 


Voi(yc) = llfj L 

\j 

= {27tr[ fir,dr,. 

The integral can now be written as 


(37) 


/t-t f'-Lj>2n /-I-E^so 

Vol(Vc) = f Hy]] (2;r)"y^ ri^ri radrs ... . 

This can be calculated by beta functions*; 

{2nf 


Voi(y,) = 


(2n+l)!^. 


Ufr 


(38) 


(39) 


This finally means that we obtain the lower bound on the region of total survivability of a linear system with no real 
eigenvalues and all real parts of the eigenvalues equal to zero of 


Volffa > , 

'The detailed calculation can be found here http://math.stackexchange.eom/a/207605. (Accessed: July 15, 2016) 


(40) 
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B Relationship to Basin Stability 

Let us further consider the relationship between basin stability and survivability. Consider the union of all attractors A in X. 
Following the terminology of,^ we split the set A into desirable attractors A+ and undesirable attractors A^. Dehne and X^ 
to be the basin of attraction of A+ and A^ respectively. The basin stability of A+ with respect to some initial region X** is then 
dehned as 


Voi(x‘'nx+) 

“ Vol(XO) 


(41) 


Note, that defining basin stability requires knowledge of the respective attractor. Efficiently evaluating it numerically 
requires a criterion to evaluate whether the system will converge to a certain attractor. 

Let us assume that the sunny region cleanly separates the set of attractors, that is, there are some attractors that are entirely 
sunny and others that are entirely shaded, but none that intersect both regions. In this case, we can establish a quantitative 
relationship between basin stability and survivability. We can choose A+ = AnX+ (i.e. desirable attractors are contained 
inside the sunny region), and we know that asymptotically the trajectories converge either to A^ or A+ or diverge. Thus every 
trajectory that does not contribute to the basin stability also has to leave the sunny region eventually and can not contribute to 
the survivability either. Thus if 


A+ =Anx+ 

and X^=X+, 

then Sb >So<,. (42) 

The difference between the two values is exactly the measure of initial conditions whose trajectories leave the sunny region 
intermittently but eventually return to it and stay. Thus we see that whether basin stability or survivability is the appropriate 
measure depends on whether the forbidden region is merely unpleasant, and we want our stay there to be hnite, or whether the 
forbidden region is deadly and we absolutely do not want the system to enter it at all. 

C Pulse-Coupled Integrate-and-Fire Oscillators 

We here give further details about the integrate-and-hre model for coupled neurons we used in the main text. Firstly, the free 
dynamics of an oscillator j is given by 


f(f) = l. (43) 

When an oscillator j reaches the threshold, 0/(f) = 1, its phase is reset to zero, = 0, and the oscillator emits a pulse 

that is sent to all oscillators i possessing an in-link from j. After a delay time T this pulse induces a phase jump in the receiving 
oscillator i according to 


t)+) :=min^l,^^^^^y-Y^-fexphe,j(j);(f4-T)^ (44) 

The phase dependence is determined by a twice continuously differentiable function t/(0) that is assumed to be strictly 
increasing, U'{^) > 0, concave (down), U''{^) < 0, and normalised such that t/(0) = 0 and U{1) = 1. 

This model, originally introduced by Mirollo and Strogatz,^ is equivalent to different well known models of interacting 
threshold elements if t/((j)) is chosen appropriately. Here we take functions of the form 

+ (45) 

where b > 0 parametrises the curvature of U, that determines the strength of the dissipation of individual oscillators. The 
function U approaches the linear, non-leaky case in the limit lim;,^ot4(0) = 0- Other nonlinear choices ofUy^Uh give results 
similar to those reported above. 

The considered graphs are strongly connected, i.e. there exists a directed path between any pair of nodes. We normalise the 
total input to each node Y!j=i £ij = £ such that the fully synchronous state exists. Furthermore for any node i all its ki incoming 
links have the same strength Eij = e/ki. 
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